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ABSTRACT 

We report observations of the transiting extrasolar planet, HD 209458 b, designed to detect 
the secondary eclipse. We employ the method of 'occultation spectroscopy', which searches in 
combined light (star and planet) for the disappearance and reappearance of weak infrared spectral 
features due to the planet as it passes behind the star and reappears. Our observations cover two 
predicted secondary eclipse events, and we obtained 1036 individual spectra of the HD 209458 
system using the SpeX instrument at the NASA IRTF in September 2001. Our spectra extend 
from 1.9 to 4.2 /xm with a resolution (A/AA) of 1500. We have searched for a continuum peak 
near 2.2 /j,m (caused by CO and H 2 absorption bands), as predicted by some models of the 
planetary atmosphere to be ~ 6 x 10~ 4 of the stellar flux, but no such peak is detected at a 
level of ~ 3 x 10 -4 of the stellar flux. Our results represent the strongest limits on the infrared 
spectrum of the planet to date and carry significant implications for understanding the planetary 
atmosphere. In particular, some models that assume the stellar irradiation is re-radiated entirely 
on the sub-stellar hemisphere predict a flux peak inconsistent with our observations. Several 
physical mechanisms can improve agreement with our observations, including the re-distribution 
of heat by global circulation, a nearly isothermal atmosphere, and/or the presence of a high cloud. 



1. INTRODUCTION 

The discovery of the first transiting extrasolar 
planet, HD 209458b, (Charbonneau et al. 2000; 
Henry et al. 2000) has led to several new observa- 
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tions designed to characterize the physical prop- 
erties of the planet. These observations have pro- 
vided a determination of the planetary and stellar 
radii and the true planetary mass (Brown et al. 
2001; Cody & Sasselov 2002). Furthermore, given 
the geometry of the orbit, we are now beginning 
to learn about the structure of the planet's at- 
mosphere. The atmosphere was first probed con- 
clusively by Charbonneau et al. (2002), who re- 
ported a detection of the sodium doublet in trans- 
mission as the planet crossed in front of the star. 
Although typical models account for the strong 
stellar irradiation and estimate that the effective 
temperature of the planet is T ~ 1100-1800 K 
(e.g., Seager & Sasselov 1998; Charbonneau et al. 
2000), no measurements of the temperature are 
available from actual observations of the planet. 
An attempt to detect the reflected starlight from 
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the planet during secondary eclipse (i.e., the time 
when the planet disappears behind the star) has 
recently been reported (Kenworthy & Hinz 2003), 
but these observations were performed in the vis- 
ible region and were limited by instrumental and 
atmospheric effects. Pioneering infrared observa- 
tions (Lucas & Roche 2002) did not achieve suffi- 
cient sensitivity to detect realistic planetary mod- 
els of HD 209458 b and similar extrasolar planet 
systems. The infrared observations reported here, 
however, have sufficient sensitivity to detect the 
thermal emission spectrum from HD 209458 b at 
the level predicted by several models for the plan- 
etary atmosphere. 

In paper I (Richardson et al. 2003), we applied 
the method of occultation spectroscopy observa- 
tions of the secondary eclipse from the Very Large 
Telescope (VLT), and using this technique over a 
narrow bandpass near 3.6 /urn, we were able to 
place limits on the abundance of methane in the 
planetary atmosphere. The most stringent limit 
applied only to an exceptionally clear model atmo- 
sphere, and only for certain values of the eclipse 
timing, which is uncertain by as much as ~ 30 min- 
utes, given the la error in the eccentricity of the 
orbit. In this paper, we report a further attempt 
to detect the secondary eclipse of HD 209458 b us- 
ing the SpeX instrument at the NASA Infrared 
Telescope Facility (IRTF). We took a different ap- 
proach for these observations; we used a broader 
wavelength range, in order to look for the infrared 
continuum of the planet near 2.2 ^m. By observ- 
ing the combined light from the star and planet, 
we search for a change in the shape of the spectrum 
as the planet disappears behind the star and later 
reemerges. The nature of this technique makes 
our observations quite sensitive to the tempera- 
ture gradient in the planetary atmosphere. 

2. OBSERVATIONS 

We obtained a total of four nights of data from 
the SpeX instrument (Rayner et al. 2003) at the 
NASA IRTF located on Mauna Kea in Hawaii. 
These are UT 19, 20, 26, and 27 September 2001, 
and secondary eclipse events were predicted for 
UT 20 and 27 September. 

The SpeX instrument is a cross-dispersed 
echelle spectrometer capable of imaging and spec- 
troscopy in the wavelength region between 0.8 and 



5.5 /Ltm with low to moderate spectral resolution 
(A/AA < 2500). We operated the instrument in 
1.9 to 4.2 fim mode using the 0.5 arcsec slit, giv- 
ing a spectral resolution of 1500; this mode gives 
nearly-continuous wavelength coverage in this re- 
gion, over six orders of the echelle. The weather 
was excellent (low water vapor), and the seeing 
conditions were good — less than 1.0 arcsec on all 
four nights. Both eclipses occurred within an hour 
of the transit of the star across the local meridian, 
meaning that the eclipse was observed with the 
object directly overhead. 

The use of a comparison star is an impor- 
tant aspect of our observational technique, and 
we selected HD 210483, a close spectral match 
to HD 209458 and nearby in the sky (see Ta- 
ble 1 in Paper I). The comparison star was ob- 
served immediately following each observation of 
the HD 209458 system. We nodded the telescope 
between the 'a' and 'b' positions on the slit to re- 
move the terrestrial atmospheric background, and 
we recorded spectra in the sequence 'abba.' A 
given spectrum at cither slit position was recorded 
with an integration time of 6 seconds per coadd 
and a total of 3 coadds; since the two objects were 
nearly the same visual magnitude, we used the 
same integration time for observing both objects. 
Typically we would record two successive 'abba' 
sets of HD 209458 and then switch to the compar- 
ison star for two more 'abba' sets, giving approx- 
imately a 50% duty cycle. It required only about 
15 seconds to switch between the two stars. Over 
the four nights of observations, we obtained 1036 
individual spectra of HD 209458 and 868 spectra 
of HD 210483, with a typical signal-to-noise ratio 
of ~ 100 near 2.2 /im. Finally, we take calibration 
spectra approximately once per hour. To account 
for flexure and other changes over the night, im- 
portant for the large SpeX instrument, we record 
flats using an IR continuum lamp, and we record 
arc lamp spectra for wavelength calibration. 

3. ANALYSIS 

In this section we discuss the details of the 
analysis process used to interpret the observa- 
tions. This process consists of rejecting outlying 
points (or 'hot pixels'), extracting the spectra from 
the raw data, and subtracting a suitable compar- 
ison spectrum to remove the telluric features. We 
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search for changes in the resulting 'difference spec- 
trum' that are synchronous with the secondary 
eclipse and are therefore due to the spectrum of 
the planet. 

3.1. Spectral Extraction 

The first step in the analysis was the rejection of 
energetic particle events and other intermittently 
bad pixels from the raw data. This was accom- 
plished using a median filter over a set of eight raw 
data frames, corresponding to two 'abba' sets. 

For extracting the spectra from the raw data 
frames, we used the Spextool software package 
(version 3.0), written by Mike Cushing and Bill 
Vacca for analysis of SpeX data (Vacca et al. 
2003). The IDL widgets allow the user to enter 
an 'ab' pair of raw images, with the corresponding 
calibration files (flats and arcs). Given an individ- 
ual 'ab' pair, the program extracts the sky back- 
ground using the region of the slit between the 'a' 
and 'b' object positions, and removes this from the 
extracted 'a' and 'b' spectra. Further, it checks for 
the residual sky background in the pair-subtracted 
image, and removes it if necessary. Currently, the 
program outputs a standard extraction of the 'a' 
and 'b' spectra; that is, it performs a sum over 
the spatial dimension to calculate the flux at each 
wavelength point. Spextool also performs a wave- 
length calibration of the extracted spectra using 
the arc frames (recorded by taking an image of 
an Argon lamp) as well as tabulated information 
on telluric lines and comparing with the extracted 
spectra. Finally, Spextool corrects for slight de- 
tector non-linearities in the SpeX array. 

3.2. Corrections to the Extracted Spectra 

The spectra of each object for a given night 
are first interpolated onto a constant and uniform 
wavelength scale. At this point we perform a sec- 
ond quality control check to find and correct any 
remaining outlying points due to uncorrected 'hot 
pixels.' For each wavelength point we perform a 
median filter over the time series of all values of 
that point for a given object during the night. A 
point is considered to be an outlier if the difference 
between the value of the point and the median 
value is greater than 4cr, and it is then replaced 
by the median value. This process effectively re- 
moves outlying points from the spectra, but it is 



only necessary to correct about 0.3-0.4% of the 
total number of points for each object. 

Next we remove from the stack any spectra that 
are clearly discrepant. The rejected spectra cor- 
respond in most cases to observations that were 
listed in our notes as questionable, either because 
of short-term changes in terrestrial atmospheric 
conditions or poor focus of the telescope. The 
stacks of spectra for both objects now contain only 
those spectra that will later be used in the calcu- 
lation of the difference spectra. 

With the individual spectra (again, of a given 
object) now interpolated onto the same wave- 
length scale and cleaned of outlying points, we 
then correct the data at each wavelength point for 
air mass. We fit a line to the log of the intensity as 
a function of air mass, as described in Richardson 
ct al. (2003, Equation 3), although in this case, we 
correct to air mass of unity. In order to ensure that 
we do not remove the effect of the planet from the 
HD 209458 spectra, we experimented with using 
the air mass correction calculated for the compar- 
ison star to correct the HD 209458 spectra. We 
found that this did not have a large effect on the 
resulting difference spectrum, showing that cor- 
recting the HD 209458 and HD 210483 spectra sep- 
arately did not introduce a bias in the resulting 
residual spectrum. An example of extracted spec- 
tra (only order 4, -2.0-2.4 fim) of HD 209458 and 
HD 210483 after the air mass correction is shown 
in Figure 1. 

We then calculate an average air mass-corrected 
spectrum of each object for that particular night. 
By comparing the average spectra of HD 209458 
and HD 210483 we can identify stellar lines; at 
this resolution (A/AA = 1500), only strong stel- 
lar lines are important. Because of the large dif- 
ference in the radial velocities of the two stars, 
|Au| ~ 55.2 km s" 1 (Nidever et al. 2002; Wilson 
1953), a given stellar line will appear at a slightly 
different wavelength in one star relative to the 
other. These show up clearly by subtracting the 
two average spectra. We remove them by linearly 
interpolating between the points on either side of 
the line. In comparing HD 209458 and HD 210483, 
we found a total of 16 stellar lines, identified as Mg 
and Si, as well as a single H line and a single Na 
line. 
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3.3. Difference Spectra 

With the groups of spectra cleaned and ad- 
justed for air mass, we are now ready to consider 
individual spectra. First, we remove the known 
stellar lines from the individual spectra in the 
same way that they were removed from the av- 
erage spectra. Note that at this point in the anal- 
ysis, the two sets of spectra (corresponding to the 
two stars) have been extracted and interpolated 
onto separate wavelength scales. Next, an opti- 
mum shift value in wavelength is calculated (for 
both the average as well as the individual spectra 
separately) by minimizing the standard deviation 
of the difference in intensity between a target and 
a normalized comparison spectrum for each order 
separately. Using the calculated wavelength shift, 
we can then interpolate each comparison spectrum 
onto the wavelength scale of each HD 209458 spec- 
trum. This ensures that the telluric absorption 
features line up in the average and individual spec- 
tra of both objects. 

Next, we calculate difference spectra by com- 
paring individual spectra of HD 209458 and HD 210483. 
As described, we record spectra of HD 209458, typ- 
ically in groups of 8, or two 'abba' sets, and then 
switch to the comparison star. The subsequent 
observations of the comparison star are recorded 
at nearly the same air mass as those of HD 209458. 
We then calculate the normalized 'difference spec- 
trum' di from 



where £, represents a single spectrum of HD 209458, 
c, is the corresponding comparison star spectrum, 
fi and Qi are normalization factors, and c is the av- 
erage comparison star spectrum for a given night 
of observations. We calculate a difference spec- 
trum for each individual spectrum of HD 209458 
by subtracting the corresponding subsequent com- 
parison spectrum; for example, the first 'a' spec- 
trum of HD 209458 in a given 'abba' set would be 
compared to the first 'a' spectrum in the subse- 
quent comparison set. The normalization factors 
are used to place the comparison spectra (average 
and individual) on the same scale as the spectra 
of HD 209458. Although the two stars are nearly 
the same brightness, we nevertheless expect that 
they will not have exactly the same intensity lev- 
els, as seen in the example in Figure 1, due to 



slit losses, variable atmospheric absorption, and 
changes throughout the night. Both normaliza- 
tion factors are calculated for the entire spectrum, 
rather than for each order independently, to ensure 
that any overall slope in the planetary spectrum 
is not removed. The factors are computed by en- 
forcing the condition that the total of all orders in 
the comparison spectrum is equal to the total of 
all orders in the spectrum of HD 209458, as in 

fi = fi^ (2) 

and similarly for gi. 

This process effectively and completely removes 
the terrestrial absorption lines that dominate the 
object spectra. Furthermore, the process also has 
two desirable side-effects. First, it removes any 
time- varying changes in the detector response that 
may not have been corrected by the flat-fielding. 
Second, it removes many of the variations in the 
continuum and line absorption that were not cor- 
rected by the fit to air mass, because these varia- 
tions may also be time- variable. Note also that the 
resulting normalized difference spectrum contains 
the candidate planetary spectrum in flux units rel- 
ative to the stellar spectrum, or the 'contrast' (Su- 
darsky et al. 2003). 

3.4. Averaging over Wavelength 

As described above, the individual difference 
spectra di have been calculated by subtracting an 
individual (i.e., at a given slit position 'a' or 'b') 
HD 210483 spectrum from the corresponding indi- 
vidual HD 209458 spectrum, as indicated by Equa- 
tion 1. We then proceed to average the individ- 
ual difference spectra di over wavelength in order 
to improve the signal-to-noise ratio. We separate 
the spectra into sections, or 'bins', and average 
the points within each bin. The bins are defined 
specifically for each region, such that a bin bound- 
ary does not fall on a telluric or stellar feature. 
This ensures that we do not introduce a bias in 
the bin average. The width of each bin is there- 
fore variable but is typically ~ 0.025 fim. 

In looking at the series of binned difference 
spectra, we see that most of them exhibit an over- 
all slope. The variation in intensity is a slowly- 
varying function of wavelength, amounting to typ- 
ically 1-5% from 2-4 ^m; a slight slope is evident 
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in the example difference spectrum from UT 27 
September shown in Figure 2. We attribute this 
effect to image motion on the slit due to guiding 
errors, coupled with the A -0 2 wavelength depen- 
dence of seeing (Linfield et al. 2001). Our inter- 
pretation is based on three considerations: 1) the 
extreme slopes are clearly related to image quality, 
as judged by the slit losses through total intensity; 
2) we constructed a simple model of the process 
and found that we expect to see such slopes, based 
on the wavelength dependence of seeing combined 
with guiding errors of plausible magnitude; and 3) 
the IRTF staff also see the effect (M. Cushing, pri- 
vate communication, 2002). We have investigated 
several methods for removing this baseline effect, 
but we found that our results for the short wave- 
length region (2-2.5 /im) are robust and do not 
change significantly based on whether a correction 
is made or on the details of the correction. This 
makes sense, given that the flux peak in the plan- 
etary spectrum predicted by most models (e.g., 
Sudarsky et al. 2003) is a relatively sharp peak 
and therefore insensitive to the removal of a low- 
order polynomial over the entire 2-4 /im region. 
We therefore present our results with no attempt 
to remove these gradual slopes. 

Finally, the uncertainties have been carefully 
propagated through the entire analysis. At the 
beginning of the analysis, after the spectra have 
been extracted from the raw frames using Spex- 
tool, we calculate the noise level in each data set. 
(Each data set consists of a time-series of spec- 
tra of a given object for a given night of obser- 
vations.) We take the error in each point to be 
the standard deviation in the time series at each 
wavelength point independently. This error value 
is propagated through the calculation of the dif- 
ference spectra and the wavelength binning using 
the standard error propagation formulae. Further- 
more, we tested our analysis technique by adding 
a synthetic planetary signal to the data and veri- 
fying that our algorithms could extract it with the 
correct amplitude. 

3.5. Fit to Eclipse Curve 

The final step in the analysis is the comparison 
of the processed stack of residual difference spec- 
tra to the eclipse timing curve. The eclipse curve 
is calculated based on the predicted time of cen- 
ter of secondary eclipse (corrected for the Earth- 



Sun light travel time), using known values of the 
period and time of center eclipse from Schultz 
et al. (2003) (P = 3.5247542 ± 0.0000044 days 
and T c = 2452223.896173 ± 0.000086 HJD) 
and assuming a circular orbit (eccentricity e = 
0). We also note that updated ephemeris data 
(P = 3.5247501 ± 0.0000004 days and T c = 
2452618.66888 ± 0.00010 HJD) is available from 
Wittenmyer et al. (2002). The eclipse curve is 
constructed in a simple way; zero represents dur- 
ing eclipse, unity represents out of eclipse, and we 
perform a linear interpolation between the ingress 
and egress points. Then we perform a linear least 
squares fit of the difference spectra to the eclipse 
curve, at each wavelength bin. The amplitude 
from the least squares fit effectively produces a 
final spectrum that represents the out-of-eclipsc 
spectra minus the in-eclipse spectra, and there- 
fore represents the candidate planetary spectrum. 

Finally, we consider the possibility of a shift in 
the predicted time of center of secondary eclipse, 
caused by a non-zero value of the orbital eccen- 
tricity (Charbonneau 2003). The predicted times 
are based on the assumption of zero eccentric- 
ity, which implies that the orbit is circular and 
that the secondary eclipses occur exactly halfway 
between primary eclipses. However, the current 
value of the eccentricity, based on a single-planet 
fit to the Doppler velocity measurements, is e = 
0.0281 ± 0.0120 with uj = 69.78 ± 1.2° (G. Laugh- 
lin, private communication, 2003). If this eccen- 
tricity is taken at face value, it implies a shift in 
the timing of the secondary eclipse of St = 31.4 
minutes, a slightly later eclipse. As shown below 
in Section 5, we find certain models for the plane- 
tary atmosphere to be inconsistent with the data 
at e = 0. We have investigated the effect of chang- 
ing the time of center eclipse by as much as 120 
minutes in either direction. At no eclipse times do 
we find evidence for the rejected model spectra at 
full strength in the data. 

The final residual spectrum shown in Figure 3 
represents the resulting fit to the eclipse curve, av- 
eraged for both of the in-eclipse nights, assuming 
zero orbital eccentricity. We show only the re- 
sult for the short wavelength region (1.9-2.5 /im). 
The larger error bars in the 3.0-4.0 /im region 
(suggested by Figure 2) are caused by the ter- 
restrial background and prevent a conclusive in- 
terpretation of the planetary spectrum in this re- 



5 



gion. Analysis of these long-wavelength data is 
continuing but is beyond the scope of this paper. 
Also plotted in Figure 3 is a model spectrum by 
Sudarsky et al. (2003) — their baseline model for 
HD 209458 b; this model is excluded by our data, 
as described in detail in Section 5. 

We note that a flat line would be consistent 
with our data shown in Figure 3, implying that 
any model that does not exhibit the peak shown 
in the Sudarsky et al. (2003) baseline model would 
escape detection in the present analysis. A flat 
line drawn through the data overlaps 16 of the 27 
points within the individual lcr error bars, with 
a reduced chi-squared value of 2.35 (suggesting 
that the error bars are slightly underestimated and 
smaller than the scatter in the data). Note that 
for a normal distribution of errors, 19 of 27 points 
would be expected to fall within ±1ct; therefore, 
our data are roughly consistent with random devi- 
ations from a flat line. We verified that the errors 
in each wavelength bin are normally distributed. 
In order to explore the range of models that are 
rejected by our data, we do some diagnostic model 
calculations. 

4. MODEL CALCULATIONS 

Our model calculations are intended to help in- 
terpret the null result for the planetary spectrum 
in the region between 2.0-2.5 /im. We have im- 
plemented a simple spectral synthesis code to cal- 
culate model spectra for HD 209458 b. The algo- 
rithm consists of adopting a temperature-pressure 
profile and computing the emergent thermal emis- 
sion spectrum. We then scale the temperature- 
pressure profile in an ad hoc fashion and re- 
compute the spectrum. Although we do not en- 
force radiative equilibrium when scaling the pro- 
file, the method is nonetheless consistent with our 
intent, which is to provide a diagnostic of the plan- 
etary atmosphere. 

The three fiducial temperature-pressure pro- 
files, described below, were used as input to the 
spectral synthesis code and were calculated in ra- 
diative equilibrium using improvements on the 
method given by Seager & Sasselov (1998) and 
Seager et al. (2000). We devised a simple way of 
scaling the fiducial temperature-pressure profiles 
to see how these scalings changed the resulting 



spectrum. The prescription we used was 

T[ =T i + 7 (T i -T )+T offsct , (3) 

where Ti is the temperature of layer i, and To is 
the boundary temperature, i.e., the temperature 
at small optical depth. Two parameters control 
the shape of the synthetic profile: 7 determines 
the temperature gradient, or atmospheric heating, 
and T ffsot is used to adjust the overall tempera- 
ture of the profile. We characterize the shape of 
the resulting profile using the temperature at op- 
tical depth unity for A = 2.15 /im (the center of 
the bandpass) and subtracting the boundary tem- 
perature: 

AT = T{t Uc = 1)-T . (4) 

Armed with a temperature-pressure profile, we 
continue by calculating the continuous opacity 
due to collision-induced absorption (CIA) for H2- 
H2 (Borysow 2002) and H2-HC (J0rgensen et al. 
2000). We include line opacities for H2O, CO, and 
CH4. The water lines are taken from the extensive 
line database calculated by Partridge & Schwcnkc 
(1997) and compiled by Kurucz, 6 the CO lines are 
from Goorvitch (1994), and the CH4 lines were ob- 
tained from HITRAN (Rothman et al. 1998). The 
relative mixing ratios of the three species were de- 
termined using the simple thermochemical equi- 
librium formulae provided by Burrows & Sharp 
(1999). We then compute emergent intensity as 
the formal solution to the radiative transfer equa- 
tion, with the source function equal to the Planck 
function, from 

^=h-B v (T). (5) 

We then solve for I v and calculate the flux density 
by integrating over /x, as in 

F v = f f B v {T(T))e- T/ »dTdii. (6) 

Finally, we calculate the contrast ratio by dividing 
by a Kurucz model (Kurucz 1992) for HD 209458. 
The parameters used in the Kurucz model were 
T = 6000 K, log 3 = 4.25 (cgs), and [Fe/H]= 0.0, 
as calculated by Cody & Sasselov (2002) using the 
parameters in Mazeh et al. (2000). 

http: //kurucz .harvard. edu/molecules/h2o/ 
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The assumption that the source function is 
given by the Planck function amounts to ignor- 
ing scattering processes in the model atmosphere. 
Only clouds and aerosols will produce significant 
scattering at this wavelength. Any attempt to 
specify the distribution of clouds in our suite of 
diagnostic models would be problematic, since the 
structure of these models has been perturbed in 
an ad hoc fashion. Fortunately, the amplitude of 
the flux peak at this wavelength is determined pri- 
marily by the wavelength distribution of CO and 
H 2 opacity and the temperature profile of the 
atmosphere. 

It is useful at this point to clarify our usage 
of the Sudarsky et al. (2003) baseline model. We 
use their tabulated flux values 7 for the planet, di- 
vided by the Kurucz model described above. We 
also multiply by the area ratio 0.016, obtained by 
taking R p = lA2Rj (Cody & Sasselov 2002) and 
i?» = 1.146 Rq. This yields approximately twice 
their tabulated contrast values, because their con- 
trast is phase-averaged (Sudarsky, private commu- 
nication, 2002). Using their temperature-pressure 
profile as input, we calculate a flux peak of similar 
but slightly larger magnitude with our simplified 
spectral synthesis code, and the differences could 
be explained by the fact that we ignore their cloud 
opacities in our calculation. Note that our result 
is based on the Sudarsky et al. (2003) baseline 
flux divided by the Kurucz model, not on our re- 
calculation of the spectrum. 

The assumed redistribution of incident stellar 
radiation is also different among the models, and 
we consider this factor in comparing with our data. 
Recalling the equation for the effective tempera- 
ture of the planet (Guillot et al. 1996), we have 

Tc^T^lHl-A)} 1 / 4 , (7) 

where T* is the effective temperature of the star, 
R* is the stellar radius, D is the orbital radius of 
the planet, and A is the Bond albedo. The pa- 
rameter / quantifies the redistribution of incident 
stellar radiation; the / = 1 case represents the sit- 
uation in which the incident radiation is evenly re- 
distributed over the entire planet, while the / = 2 
case represents the one in which the incident radia- 
tion is completely absorbed and re-emitted on the 

http : //zenith . as . arizona . edu/~burrows/ 



day side of the planet. The Sudarsky et al. (2003) 
baseline model assumes / = 2 (or equivalently, 
/ = 0.5 using their definition of / (Saumon et al. 
2002)). As noted, our first case for comparing with 
the observational results is the baseline Sudarsky 
et al. (2003) flux divided by the Kurucz model. 
We choose three other cases as fiducial models for 
comparison with the data. These have been se- 
lected in order to test a wide range of the pos- 
sible physical conditions in the atmosphere. We 
choose two / = 1 models based on the formula- 
tion by Seager & Sasselov (1998); Seager et al. 
(2000), one with no clouds and one with cloud 
opacity due to MgSi0 3 , Fc, and A1 2 3 . Since the 
Sudarsky et al. (2003) baseline model considers 
cloud opacities and assumes / = 2, we further se- 
lect a cloudless model with / = 2 (based on the 
Seager & Sasselov (1998); Seager et al. (2000) for- 
mulation). The temperature profiles of these four 
cases are shown in Figure 4, and the results for the 
four flux peaks (the Sudarsky et al. (2003) baseline 
result as well as the cloudless and cloudy fiducial 
cases) are shown in Figure 5. 

5. RESULTS AND INTERPRETATION 

We now fit these four modeled flux peaks to the 
data using linear least squares, as in 

Ti = a + brrii, (8) 

where represents the resulting difference spec- 
trum (e.g., Figure 3), m t represents the contrast 
ratio from a given model, i represents the wave- 
length bin, and a and b are the coefficients of the 
fit. The slope of this fit, b, which wc call the 'model 
amplitude,' therefore represents the 'amount' of 
the model signal that appears in the data. A 
model amplitude of unity would mean that the 
model clearly fits the data (depending on the er- 
rors), while a value of zero represents no correla- 
tion between the model and the data. Note that 
the y-offset between the model and the data is 
irrelevant, since we are searching only for the ev- 
idence of the shape of a given model within the 
data. The results for the four fiducial models are 
shown in Table 1. For each model, we indicate 
the model amplitude 6, the uncertainty in b, and 
the reduced chi-squared of the fit. We further list 
the reduced chi-squared for the case in which we 
force the model amplitude to be unity; that is, we 
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assume that the model appears in the data at full 
strength. For all four fiducial models, the slope 
of the fit is roughly zero, suggesting that none of 
them fit the data. However, the rejection is much 
stronger for the both / = 2 models, as seen by the 
reduced chi-squared values for forcing the model 
amplitude to be unity. These high values indicate 
that it is extremely unlikely that the / = 2 mod- 
els fit the data. Given that the observations and 
analysis are highly differential, we are furthermore 
confident that a peak in the planetary spectrum is 
not being fortuitously canceled by some system- 
atic error in the data. 

The least squares fit described above is diffi- 
cult to interpret for models with small amplitude 
flux peaks, as in some of our scaled models (see 
Equations 3 and 4). ft is also desirable to explore 
another way to quantify the degree to which the 
four specific models are rejected by the data. We 
use a Monte Carlo technique to accomplish this. 
To quantify the statistical significance of the rejec- 
tion of any particular model, we constructed 'fake 
data' by adding normally-distributed noise to the 
modeled contrast values, with o equal to the size 
of the error bar in each bin of Figure 3. We cre- 
ated f 5 fake data sets for each model and then fit 
each one to the original model using linear least 
squares (in the same way in which the real data 
were fit to the model to calculate the model ampli- 
tude). Next, we determine the number of resulting 
model amplitudes which are greater than the value 
obtained by fitting the original model to the actual 
data; this is a liberal criterion for concluding that 
the model is consistent with a given fake data set. 
The percentage of fake data sets consistent with 
a given model, using this criterion, is the 'detec- 
tion efficiency', Q. A high value of Q, say > 99%, 
implies that a given model consistently produced 
at least a minimal signature in the fake data. The 
results for the four specific models discussed thus 
far are shown in Table 2, and are also plotted in 
Figure 6 (filled symbols). A model with a high de- 
tection efficiency using the Monte Carlo technique 
implies that it is detectable. Thus, the two / = 2 
models — having model amplitudes near zero and 
high detection efficiencies — are again concluded to 
be rejected. These results are consistent with the 
chi-squared analysis described above. 

We tested our suite of scaled models, with the 
prescription given by Equation 3, using this tech- 



nique. The parameters and results of these mod- 
els, as well as the model amplitudes, are shown in 
Table 3 for the cloudless / = 1 cases and in Table 4 
for the cloudy / = 1 cases. The results of Q vs. 
AT are plotted in Figure 6 for both sets of models. 
We reject a model if it has a value of Q > 99% and 
also has a model amplitude near zero. Looking at 
Tables 3 and 4, the rejected models correspond to 
the steeper T/P profiles, with 7 > 1.25. Thus, we 
can sec that the temperature profiles with a higher 
temperature gradient, causing a larger peak near 
2.2 /im in the emergent spectrum, can be elimi- 
nated as inconsistent with our observations. The 
more isothermal profiles, on the other hand, are 
consistent with the data. Figure 6 suggests that 
only models with AT < 1100 K arc consistent 
with our data. Although the rejection limit also 
depends on T g sct from Equation 3, it is strongly 
dependent on AT. 

6. CONCLUSION 

Our results are extremely sensitive to the 
temperature-pressure profile, which is as expected, 
given that the Planck function varies strongly with 
temperature in this wavelength region. Besides 
the shape of the profile, there are a variety of 
other possible processes in the planetary atmo- 
sphere that could render the planet undetectable 
with these observations. For example, the broad 
peak in the spectrum near 2.2 /xm is caused by the 
large number of weak water lines on either side of 
this feature. Absorption by these lines lowers the 
continuum level on either side, thus creating the 
peak. One obvious process by which the peak 
would not appear is a lowering of the water abun- 
dance in the planetary atmosphere. We tested this 
idea and found that only when the water abun- 
dance was lowered by a factor of 10 did the peak 
become small enough to escape detection using 
our method. This value for the water abundance 
(~ 9.8 x 10~ 5 ) seems unrealistically small. 

We suggest three ways in which the 2.2 /jm flux 
peak would be lowered sufficiently to be consistent 
with our observations. The first is the efficient re- 
distribution of the incident radiation over the en- 
tire planet, which can be achieved through atmo- 
spheric dynamics, i.e., winds. Depending on how 
efficiently the incident radiation is redistributed, 
the effect could be a temperature asymmetry be- 
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tween the 'day' and 'night' sides of the planet. 
Showman & Guillot (2002) suggest that the day- 
night temperature difference may be as much as 
300 K; even a difference of this magnitude would 
lead to strong winds of 1 km s _1 or more. The 
/ = 2 situation would tend to steepen the T/P 
profile in the atmosphere, since the strong stellar 
irradiation would have to be absorbed relatively 
deep in the atmosphere and re-radiated. The two 
specific f = 2 models considered above, which as- 
sume that the heat is absorbed and re-radiated 
only on the day side, are inconsistent with our 
data (as shown in Table 1). This suggests that 
the large temperature asymmetry scenario may 
not be viable and that some sort of transport 
mechanism may be at work to move the incident 
energy from the hot to the cold side. We cau- 
tion that we have observed only two secondary 
eclipses of HD 209458 b; recent simulations by Cho 
ct al. (2003) indicate that the temperature asym- 
metry between the two hemispheres may be time- 
variable. This effect can be so strong as to render 
the day side to be colder than the night side at 
times. 

A second mechanism by which the 2.2 /im flux 
peak could be lowered is additional opacity in the 
upper atmosphere, which would flatten the T/P 
profile and thus reduce the relative strength of 
the peak. For example, the T/P profiles shown 
by Barman et al. (2001) are nearly isothermal, 
and will certainly produce 2 /zm spectra consistent 
with our observed limit. Barman (2003, private 
communication) attributes this behavior to their 
treatment of the opacity, specifically the inclusion 
of strong gaseous opacity due to metals, which dif- 
fers from the opacity treatment in other models. 
However, consistency with our data does not nec- 
essarily indicate that a given model is correct, and 
we leave it to the theorists to sort out the differ- 
ences among competing models. We suggest that 
a high-altitude cloud (first suggested by Seager 
& Sasselov (2000) and favored by Charbonneau 
et al. (2002) to explain their low sodium result) 
would also make the atmosphere more isothermal 
if the cloud is sufficiently absorbing and vertically 
extended. Note, however, that the Charbonneau 
et al. (2002) sodium result requires placing such 
a cloud at the limb of the planet, whereas our re- 
sults require the clouds to be broadly distributed 
on the day side. 



A high cloud might also affect the emergent 
spectrum through a third mechanism. For a suffi- 
ciently high cloud, the atmosphere above the cloud 
deck would be thin, giving a small column den- 
sity in which a spectral signature could be formed. 
Very little incident radiation would penetrate an 
optically thick cloud, meaning that the tempera- 
ture structure would be dominated by the black- 
body reflection or emission of the cloud itself. The 
emergent spectrum would be determined by the 
small part of the atmosphere above the cloud deck, 
rather than the details of the T/P profile. 

Finally, we note that our results have been ob- 
tained using the 3-meter 1RTF, indicating that 
valuable observations of extrasolar planets, lead- 
ing to new information on the atmospheres of 
these objects, can be conducted using modest- 
aperture telescopes from the ground, given favor- 
able observing conditions and careful analysis. 
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is supported in part by a NASA Graduate Stu- 
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the Office of Space Science at NASA Headquar- 
ters (grant number NGT5-50273). Other aspects 
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gins of Solar Systems Program. Sara Seager is 
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Table 1 

Results of the least squares fit of the data to each of the four fiducial models, 
indicating the slope (model amplitude) b and the uncertainty in b, as well as the reduced 
chi-squared for the fit. the rightmost column indicates the reduced chi-squared value 
for the case in which we force the slope to be unity. 



Model 


b 




x 2 


X 2 (b = 1) 


Cloudless (f=l) 


0.02 


± 0.29 


2.35 


2.85 


Cloudless (f=2) 


0.03 


± 0.10 


2.34 


6.43 


Cloudy (f=l) 


-0.25 


± 0.24 


2.30 


3.43 


Cloudy (f=2)* 


-0.01 


± 0.14 


2.35 


4.66 



*As calculated by dividing the flux from the Sudarsky 
ct al. (2003) baseline model by the Kurucz model for the 
stellar flux; sec Figure 5 (thin, middle curve). 
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Table 2 

Results from the Monte Carlo technique for the four fiducial models. The value of AT 
from Equation 4 is shown, as well as the value Q describing the efficiency of the 

DETECTION USING THE MONTE CARLO TECHNIQUE, AS PLOTTED IN FIGURE 6. 



Model 


AT (K) 


Q (%) 


Cloudless (f=l) 


960.36 


99.3 


Cloudless (f=2) 


1265.38 


100.0 


Cloudy (f=l) 


797.92 


100.0 


Cloudy (f=2)* 


~ 1180 


100.0 



* As calculated by dividing the flux from 
the Sudarsky ct al. (2003) baseline model 
by the Kurucz model for the stellar flux; 
sec Figure 5 (thin, middle curve). 
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Table 3 



Results for scaled cloudless / = 1 profiles. The first column is the value of AT from 
Equation 4, and the scaling parameters are defined in Equation 3. The model amplitude b 
and uncertainty c7 b are also shown. the value q represents the detection efficiency 
using the Monte Carlo technique, as plotted in Figure 6 (diamonds). 



AT (K) 
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^offset 
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Q (%) 
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05 
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00 
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00 
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99.9 


1328.16 
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13 


100.0 


1300.31 
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100 





03 





11 


100.0 
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14 


100.0 
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00 
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100.0 
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00 
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100.0 
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00 


-100 





04 





08 


100.0 
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Table 4 



Results for the scaled cloudy / = 1 profiles. The same quantites as described in Table 3 
are shown. The values for the detection efficiency Q are also plotted in Figure 6 

(squares). 



AT (K) I 7 T offsct b 7 b Q (%) 



418.63 
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Fig. 1. — Sample extracted spectra of HD 209458 
(lower, thick curve) and HD 210483 (upper, thin 
curve) from UT 27 September (near air mass of 
1.27). The comparison star has not yet been nor- 
malized for subtraction from the HD 209458 spec- 
trum. The overall shape of the spectra is due to 
the instrument response function. Spectral fea- 
tures are primarily telluric. 
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Fig. 2. — Sample binned residual spectrum from 
UT 27 September, corresponding to spectra shown 
in Figure 1. Five different orders from the spectro- 
graph (three in the short wavelength region (1.9- 
2.5 /j,m) and two in the long wavelength region 
(3.0-4.1 /xm)) are shown as alternating filled and 
empty symbols. Because of large uncertainties re- 
sulting from the terrestrial background, the long 
wavelength data are not considered further in this 
analysis. 
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Fig. 3. — Final difference spectrum, calculated 
by averaging the data from the two 'in-cclipsc' 
nights. This was obtained using 550 individual 
spectra of HD 209458, with an equal number of 
HD 210483 spectra. The Sudarsky et al. (2003) 
baseline model is over-plotted at a higher spec- 
tral resolution (A/AA ~ 2200) than the binned 
data. The offset between the model and the data is 
not important, since we are comparing the shapes 
of the two spectra; thus, the mean has been sub- 
tracted from the data and model, respectively, for 
plotting purposes. 
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Fig. 4. — Temperature-pressure profiles of the four 
fiducial cases. The dashed curve is from the Su- 
darsky et al. (2003) baseline model, the thick solid 
curve is the cloudy / = 1 case, the thin lower curve 
is the cloudless / = 1 case, and the dotted line is 
the cloudless / = 2 case. 
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Fig. 5. — Comparison of contrast for the four fidu- 
cial cases. Uppermost dashed curve is the cloud- 
less / = 2 case. The thick middle curve is the 
cloudy / = 1 case, and the lower dotted curve is 
the cloudless / = 1 case; these three are the re- 
sult of the simple spectral synthesis calculations, 
plotted at a spectral resolution of 1500. The thin, 
solid curve is the Sudarsky et al. (2003) baseline 
result for the planetary flux divided by the Kurucz 
model (see text). 
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Fig. 6. — Detection efficiency Q for a suite of 
scaled models vs. AT, which represents the dif- 
ference between the temperature at optical depth 
unity at the center of the bandpass and the bound- 
ary temperature (see Equation 4). See Tables 3 
and 4 for details of each scaled profile. The filled 
symbols, along with the asterisk, represent the 
four fiducial cases described in Table 2. 



